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Abstract 

The least-squares transformation of a discrete-time multivariable linear 
system into a desired one by convolving the first with a polynomial system 
yields optimal polynomial solutions to the problems of system compensation, 
inversion and approximation. The polynomial coefficients are obtained from 
the solution to a so-called normal linear matrix equation, whose coefficients 
are shown to be the weighting patterns of certain linear systems. These, in 
turn, can be used in the recursive solution of the normal equation. 
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1 . Introduction 


The transformation of a given linear system into one having desired 
properties, by cascading it with another linear system, has been known in 
control system design as cascade compensation (see, e.g., [1]). The problems 
of system inversion and approximation can also be formulated as system trans- 
formation problems by properly defining the roles of the systems involved. An 
exact causal transformation from the given to the desired system may not exist 
or may not be implementable because of complexity constraints and an approxi- 
mate transformation of a relatively low complexity may be desired. However, 
finding an approximation within the class of rational linear systems of a 
given order is, typically, a difficult parameter optimization problem involv- 
ing local extrema (see, e.g., [2]). 

In this paper we consider the transformation of discrete time, multivari- 
able linear systems by convolving them with systems having finite weighting 
patterns or, equivalently, polynomial transfer matrices. In contrast to the 
continuous-time case, discrete time polynomial systems do not present a reali- 
zation problem. Such systems have been used extensively in statistics and in 
signal processing for predicting and filtering discrete time processes (see, 
e.g., [31, [4] ) . The proposed approach may be viewed as the compensator 
analog of Levinson's (polynomial) approximation to Wiener's (rational) filter 
[5]. The least-squares approximation criterion yields a globally optimal 
solution to the system transformation problem at hand in the form of the 
solution to a so-called normal linear matrix equation, similar to the one 
arising in the filtering problem. 

Representing the given and the desired systems in state-space, we obtain 
the coefficients of the normal equation and of the approximation error in 
closed explicit forms and show that they are the weighting patterns of certain 


linear systems. These systems can, in turn, be used in mechanizing the recur- 
sive solution of the polynomial coefficients by a multivariable version of the 
Levinson procedure [6]. The applicability of the polynomial transformation 
method to the problems of system compensation, inversion and approximation is 
discussed and illustrated by numerical examples. 

2. Least-Squares Polynomial Compensation. Inversion and Approximation 
of Linear Systems 
System Compensation 

Given two discrete-time linear systems, having rational transfer matrices 
of the same dimensions Q(z) and n(z), it is desired to find a system having a 
polynomial transfer matrix 

*(z) * 2 Ti 2 ” 1 (2.1) 

i=0 1 

that, cascaded with Q(z) will approximate n(z) in some sense. Denoting by 
{p k J, {q k J and { t^} , k 2 0, the weighting patterns corresponding, respec- 
tively, to Q(z), n(z) and *(z), an approximation measure is defined by the Lg 
norm 


IQ(z)e(z) - n(z)i 2 = |tr jf 2n [o(e J ")*(e J “) - n(e>)]* 

x [a(e > )*(e Jw ) - n(e J ")]dwJ 1/2 

* { tr £ 0 (p k * \ - \) T (p k * - v }' /2 <2 - 2) 
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where tr denotes trace, (•)* denotes complex conjugate transpose, (») T 
denotes real transpose, and * denotes the convolution operation, that is 


\ ' e k * £ 0 Pk-i'i 


Denoting 


u . rmT mT t T,T 
t » (I 0 T, ••• T n J 


it can be verified that the approximation measure (2.2) takes the form 


lo(z)«(z) - n(z)l 2 = (tr(Lt - g) T (Lt - g )} 1/2 


Pn Pn-1 
Pn + 1 p n 


9 s hi 


Minimizing (2.3) with respect to t yields the normal equation 


L T Lt = L T q 
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Denoting 



(2.5) 


and 



( 2 . 6 ) 


where 


• E 

k=0 


p k+i p k 


and 


-i 


= E 

k=0 


p k°*k-t-i 


(2.7) 


the normal equation becomes 


Rt = s 


( 2 . 8 ) 


The optimal approximation error is obtained by substituting (2.4) into 
(2.3) as 



(2.9) 
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where t^ is defined by the solution of (2.8) and 

v = f) qja (2.10) 

u k=0 K K 

Suppose that both Q(z) and n(z) are stable, then 

tr E pip,, < - and tr £ q£q k < - 

k=0 ksO K 

implying that both r^ and s^ have finite limit values and, consequently, the 
normal equation (2.8) and the optimal error (2.9) are well defined. We note 
that as the polynomial degree n is increased to infinity, the optimal error 
vanishes, provided that an exact causal linear transformation between 8( z) 
and n(z) exists. This can be seen immediately from the fact that the error 
corresponding to the optimal sequence {t^} is minimal, while the one corre- 
sponding to the impulse response of the exact transformation converges to 
zero. Suppose that for some l > 0 we have p^sO, 0 < k < i - 1 and 
p t i 0 and for some m>0, q k sO, 0<ksm- 1 and p m i 0. Then the 
systems Q(z) and n(z) are said to have, respectively, l and m delays or 
"zeros at infinity” (see, e.g., [7], p. 449). If t > m, an exact causal 
rational transformation from a(z) to n(z) clearly does not exist and the 
error corresponding to a polynomial compensator will not vanish even when its 
degree is increased to infinity. If, on the other hand, i < m, the compensa- 
tion problem can be simplified somewhat by replacing 0(z) by z“ v ; Q(z) and 
then multiplying the resulting compensator by This is illustrated 

by an example in section 4. 
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System Inversion 

The problem of system inversion is a special case of the compensation 
problem. Suppose that it is desired to find a polynomial system of the form 
(2.1) that approximates the inverse of a given rational system Q(z). An 
inversion criterion is defined by minimizing the norm 

»Q(z)*(z) - Il 2 = |tr p (p k » t k - I« k ) T (P k * t k - I« k )} 1/2 (2.11) 

which yields the normal equation 


Rt = b (2.12) 

where R is defined by (2.6) and b is the n +• 1 block matrix 

i 

b = (p Q 0 ••• 0] T 

The optimal inversion error is obtained from (2.9) as 

e = (tr( I - pJt Q )} 1/2 (2.13) 

where t Q is obtained from the solution of (2.12). In analogy with the 
equivalent problem of solving an overdetermined linear matrix equation using 
the least-squares criterion, the term pseudo- inverse seems appropriate for the 
resulting approximate inverse system. If for some i > 0 we have p k s 0, 

0 < k < t - 1 and p t i 0, the solution of (2.12) is trivial, t = 0. How- 
ever, if the purpose of the inverse system is to reconstruct the input 
sequence from the output sequence, a more sensible approach to the inversion 
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problem is to invert a shifted version of the system, having no delays. This 
means that Q(z) is replaced by z l a(z) and the sequence Cp k J is replaced in 
equations (2. 1 1 )-(2. 13) by the sequence (p£ l) ), ^ > 0, where p£ l) = P k+1 - 
The resulting approximate inverse system should then be viewed as operating on 
a record of the output sequence, advanced i time units (u k+1 replacing u k ) 
to produce the input sequence. 

System Approximation 

Suppose that it is desired to approximate a given rational system Q(z) 
by a polynomial system *(z). Writing the approximation measure as 

l«(z) - a(«)lj » {tr i t k * - P k } 

the problem can be seen to be a special case of the polynomial transformation 
problem with I and 0(z) replacing Q(z) and n(z), respectively. The result- 
ing normal equation (2.8) has 

mm T T 

R a I and s = [p Q p 1 • • • p R ] 

yielding 

t = ( P J P^ ••• pJ) T <2. I'D 

It follows that the best polynomial approximation of a linear system has the 
first markov parameters of the system as polynomial coefficients. 

When the given system has a weighting pattern of long duration (that 
is, of relatively large values for a relatively long time) a polynomial 
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approximation, which has a finite weighting pattern, will require many poly- 
nomial coefficients for a sensible representation. An alternative approach 
that does not enforce a finite weighting pattern but still maintains the 
mathematical simplicity of the polynomial transformation, is to use the inver- 
sion criterion as a system approximation criterion. This means that the 
approximate system is the one whose inverse, when convolved with the given 
system, approximates the identity matrix. The solution of the least squares 
polynomial inversion problem yields an inverse system *(z) of the form (2.1). 
The question of interest is whether *(z) has a causal inverse; that is, 
whether there exists a proper rational matrix H(z) that satisfies 

H(z)*(z) * I (2.15) 

It follows from a fundamental result for polynomial matrix equations (see, 
e.g., [7], p. 387) that a rational matrix H(z) that satisfies (2.15) exists 
if and only if Tq has full column rank and that one such matrix is given by 

H(z) = e -1 (z) (2.16) 

where e(z) consists of the elements of *(z) corresponding to a non-singular 
minor of Tq having a maximal dimension. Since e(z) is a polynomial matrix, 
its inverse is obviously proper. We note that in the single input single 
output case the approximate system (2.16) is of the autoregressive (all pole) 
type. 
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3. The Normal Equation Coefficients and its Recursive Solution 

Let (A, B, C, D) and (F, G, H, E) be state-space realizations of 
and n(z), respectively. Then, noting that 



it can be verified that 


r . 

1 


Vd * b t pb 




^B T (A T ) i_1 C T D ♦ B T (A T ) i PB 


i = 0 
i > 1 


and 


T T 

Da -*• B 1 QG i = 0 

T i-1 T i 
^THF* 'G ♦ B i QF i G i > 1 


where 


P = f) (A T ) k C T CA k 
k=0 


and 


Q = JJ (A T ) k C T HF k 
k=0 


0(z) 


(3.D 


(3.2) 


(3.3) 


(3.4) 
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the term v Q in the optimal error, defined by (2.10) can be written as 


v Q z E T E + G T VG 


(3.5) 


where 


II: j (F T ) k H T HF k 
k=0 


(3.6) 


It can be seen that the matrices P, Q and V, defined by equations (3.3) , 
(3.4) and (3.5) are central to the normal equation and the optimal error. 
Suppose that both A and F have their eigenvalues strictly inside the unit 
circle. Then it can be verified that the matrices P, Q and V are the limit 
values of the matrices P m , Q m and V m that satisfy the recursive equations 


P . s A T P A ♦ C T C 
m+1 m 


(3.7) 


v, ■ # V * cT « 


(3.8) 


V . * F T V F + H T H 
m+1 m 


(3.9) 


yielding the algebraic equations 


P s A T PA + C T C 
Q * A T QF + C T H 

and 

V = f t vf + h t h 


( 3 . 10 ) 

(3.11) 

( 3 . 12 ) 
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Equations (3.10) and (3.12) are Lyapunov equations that can be solved by 
standard methods. Equation (3.11) is not a Lyapunov equation, as the matrices 
involved are generally not symmetric. It can be solved by equating its two 
sides term by term, which results in a linear equation of the form LX = z, 
where L has dimension n^ x n^. The three equations (3- 10)— (3- 12) can be 
put in the form of a single Lyapunov equation 


P 







(3.13) 


and solved by standard procedures. These equations can also be solved by 
propagating equations (3. 7) -(3. 9) until convergence is observed. Other alter- 
natives for the solution of Lyapunov equations can be found in the literature 
(see, e.g., [8], p. 67). We note that, as can be seen from (3.3), (3.4) and 
(3.6), P has full rank if and only if (A, C) is observable, V has full rank 
if and only if (F, H) is observable, and Q has full rank if and only if both 
are observable. 

It can be seen from (3.1) and (3.2) that the coefficient sequence {r^}, 
i > 0 is the weighting pattern of the system 


Vi ■ A \ + (cT ° ♦ ftTpB)u k 

y k = B \ + (dt ° + eTPB > u k 


(3.14) 


and, similarly, the coefficient sequence {s^}, i > 0 is the weighting pattern 
of the system (F, G, D^H + B^QF,. D^E + B^QG) . 
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The symmetric Toeplitz structure of the normal equation (2.8) can be 
exploited in its solution. Employing a multivariable version of the Levinson 
algorithm (see, e.g., [6], pp. 241-244) the solution for a polynomial of 
degree n can be obtained from the solution for a polynomial of degree 
n-1. Such procedure can be used for calculating a polynomial of a given 
degree without resorting to matrix inversion, or for recursively increasing 
the polynomial order, until an error criterion is met. We note that the key 
step of the Levinson procedure is an operation of the type 


n 

£ P n+1-i fc n,i 


(3.15) 


where t_ . is the i'th parameter of a polynomial of order n and 
n f x 

r n+1-i is def>ined by (3.1). In the filtering application of the Levinson 
algorithm, both the memory requirement of the {r^} sequence and the execution 

t 

of (3.15) make the algorithm impractical for large n values. On the other 
hand, for the present application, both problems can be circumvented by noting 
that o n is the final output of the system (3.14) when the input is the 
sequence {t R q, t n ... t R n , 0}. It follows that (3.15) can be executed 
for each n by passing the polynomial coefficients of the preceding iteration 
through the system (3*14). 

4. Examples 

The numerical implementation of the proposed solutions to the compensa- 
tion, inversion and approximation problems has been examined by solving numer- 
ical examples. It should be emphasized that, in general, the polynomial 
degree needed for a satisfactory accuracy depends on the specific systems 
involved. 
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Example 1 : System Compensation 


Suppose that It is desired to compensate a system whose transfer function 
is given by 


q{2) . (2 ± 0.2)(z - 0.5 ±0.2j)(z + 0.02 ±0.1,1) 

vz; ' (z - 0.2)(z 0.1)(z + 0.6 ±0.3J)(2 - 0.1 ±0 . 3 J ) 

by a polynomial compensator , so as to approximate the desired system 

D(2) s (2 + 0.3521 )(z - 0.8521) 


The impulse responses of both systems are shown in Figure 1. State space 
representations for Q(z) and H(z) have been obtained in controller canonical 
form and the calculations described in sections 2 and 3 were carried out using 
standard software. We note that since n(z) has two zeros at infinity, it was 
multiplied by z 2 and the resulting compensator was multiplied by z" 2 , in 
order to reduce the number of optimized polynomial coefficients. The result- 
ing polynomial compensator coefficients for degrees 4, 9, 14 and 19 are listed 
below. 


II 

■=r 

4JI 

[0 

0 

0.8082 

1.6574 1.8709 

-9 3 

[0 

0 

1.0015 

2.3625 3.3676 


1.9271 

1.0806 

0.3568] T 


1.3566 0 . 5323 ] T 

3.8086 3.7758 • 3-3627 2.7198 
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-14 

= [0 0 

1.0001 

2.3603 

3.3723 

3.8409 

3.8777 

3.6051 

3.2113 


2.7782 

2.3621 

1.9735 

1.6050 

1.2353 

0.8560 

0.4754 

0. 15653 T 

and 


• 







-19 

o 

o 

l—l 

II 

1.0000 

2.3600 

3.3716 

3.8394 

3.8751 

3.6015 

3.2078 


2.7785 

2.3753 

2.0179 

1.7116 

1.4520 

1.2317 

1.0417 

0.8719 


0.7120 

0.5501 

0.3824 

0.2128 

0.0701 ] T 





The approximation errors for the four compensators were obtained from (2.9) as 
« 4 * 1.1985 , e g * 0.5726 , * l4 * 0.2544 and e ig = 0.1141 

The impulse responses of the compensated system are shown in Figure 2, 
along with the impulse response of the desired system. It can be seen that 
the compensation quality improves as the compensator's degree is increased and 
that the compensator emphasizes matching the large initial response values, as 
should be expected. 
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Figure 2 The impulse responses of the compensated system a(z)*(z) for poly- 
nomial compensators of degrees (a) 4, (b) 9, (c) 14, and (d) 19. 
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Example 2 : System Inversion and Approximation 

Inverse polynomials of degrees 3 and 6 for the system Q(z) of example 1 
were derived using the proposed method. These are given by 

* 3 (z) s 0.8493 + 1.3887Z" 1 + 1.0356z“ 2 + 0.3555z" 3 


and 

* 6 (z) = 0.9970 + 1.8479Z* 1 + 1.8620z -2 ♦ 1.3874z‘ 3 + 0.8499z" 4 + 0.3834z" 5 
♦ 0.1067z" 6 

The inversion errors associated with these and with polynomials of degrees 7 
and 9 were obtained by (2.13) as 

e 3 = 0.3888 s 0.0548 e ? = 0.0200 e g = 0.0011 

The Impulse responses of * 3 ^(z) and *^(z) are shown in Figure 3 along with 
the impulse response of Q(z). It can be seen that a close approximation is 
obtained even for these relatively low degree polynomials. A graphically 
perfect match was obtained for polynomials of degrees 7 or higher. In con- 
trast to this approximation by inversion, approximation by (2.14) can be seen 
to be inadequate. 
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Figure 3 The impulse responses of the inverse approximate (z) of Q(z) 
for polynomial degrees (a) 3 and (b) 6. 

5. Conclusion 

The least-squares parameters of a polynomial system that approximately 
transforms a discrete-time linear system into another satisfy a linear equa- 
tion, whose coefficients have been derived in terms of the state space parai* 
eters of the systems involved. The resulting method has been shown to apply 
to the problems of system compensation, inversion and approximation. 
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